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Abstract 



Mh| We present the program BOKASUN for fast and precise evaluation of the Master 

r~| I Integrals of the two-loop self-mass sunrise diagram for arbitrary values of the inter- 

nal masses and the external four-momentum. We use a combination of two methods: 
a Bernoulli accelerated series expansion and a Runge-Kutta numerical solution of a 
^ I system of linear differential equations. 
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Nature of problem: 

Any integral arising in the evaluation of the two-loop sunrise Feynman diagram can 
be expressed in terms of a given set of Master Integrals, which should be calculated 
numerically. The program provides with a fast and precise evaluation method of the 
Master Integrals for arbitrary (but not vanishing) masses and arbitrary value of the 
external momentum. 

Solution method: 

The integrals depend on three internal masses and the external momentum squared 
p^. The method is a combination of an accelerated expansion in 1/p^ in its (pretty 
large!) region of fast convergence and of a Runge-Kutta numerical solution of a 
system of linear differential equations. 

Running time: 

To obtain 4 Master Integrals on PC with 2 GHz processor it takes 3 /xs for series 
expansion with calculated in advance coefficients, 80 /US for series expansion without 
calculated in advance coefficients, from few seconds up to few minutes for Runge- 
Kutta method (depending on the required accuracy and the values of the physical 
parameters). 



LONG WRITE-UP 



1 Introduction 



The sunrise diagram with arbitrary masses is one of the basic ingredients 
of any two-loop calculation, and its fast numerical evaluation is of direct 
interest in Monte Carlo simulation programs. Many procedures for a pre- 
cise evaluation of the Master Integrals (Mi's) can be found in the literature 
[1,2,3,4,5,6,7,8,9,10,11,12,13,14]. In [15] a fast and precise series expansion was 
proposed in the much simpler case of equal internal masses. In this paper we 
adapt the method to the arbitrary mass case by considering the accelerated 
expansion in inverse powers oi p^, which provides with a fast and precise con- 
vergence in a wide region covering the biggest part of p^ values. The expansion 
does not work in a relatively small region (roughly from p^ = to the physical 
threshold). A similar expansion in p"^ around the regular point p^ = (the first 
terms where given in [16]) is unpractical due to the severe numerical insta- 
bility of the coefficients of the expansion [17], associated with the presence of 
nearby pseudothresholds (this feature is peculiar to the arbitrary mass case, 
as opposed to the equal mass limit). In the regions were the expansion in 1/p^ 
does not work we use the Runge-Kutta algorithm developed in [11,18,19] to 
obtain the Master Integrals (Mi's) of the sunrise diagram as the numerical 
solutions of a suitable system of linear differential equations. The execution 



time is of the order of a few seconds (minutes) when the Runge-Kutta method 
is used (depending on the values of p^ and the masses and also on the required 
precision), while it drops to about 80 /iS when the new expansion applies and 
to 3 /is if the expansion coefficients are calculated in advance (all CPU times 
are given for a 2 GHz PC). 



2 The notation 



We use here the same notation and definitions as in [1 1] , which we recall shortly 
for convenience of the reader. The four Master Integrals (Mi's) related to the 
general massive 2-loop sunrise self-mass diagram in n continuous dimensions 
and with fully Euclidean variables are defined as 
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where j = 0, 1, 2, 3 refers to the 4 Mi's; for j = 0, ai{j = 0) = 1 for i = 1, 2, 3; 
for j > 0, ai{j) = 2 when i = j and ai{j) = 1 when i ^ j. 

The mass scale is chosen as 

H = mi+m2 + ms , (2) 

which comes out to be the appropriate mass scale parameter for the numerical 
discussion. The expansion of the Mi's around n = 4 has the form [16] 

(3) 
where the coefficient C{n) 

C(„) = {2V^f-"'T{3-'l), (4) 



not to be expanded, can be replaced by its value C(4) = 1, at n = 4, when 
multiplying a function regular in (n — 4). The coefficients of the poles in 
{n — 4) of Fj{n, ml, m\, m^^p^) are known in closed analytic form [16,11], and 
are not reconsidered here, as from now on we deal only with the finite parts 
Ff\ml,ml,ml,p^) of the Mi's. 

It is convenient to use reduced masses and reduced external invariant 
mi + 7712 + m3 (mi + m2 + mi^y 



together with a dimensionless version of FQ{n,m\,m\,m\,p'^), defined by 

Fo^r{n,m^,m^,mr^,p) = — ■ ■ -— ; 6 

(mi +m2 + m3)^ 



as the other Master Integrals are already dimensionless, the values of all the 
functions are also dimensionless. In terms of the new variables p'^,m,i^j. the 
threshold is located at pff^^ = — 1. 

Moreover we do not write anymore, for short, the arguments of the func- 
tions and the superscript (0), so we set Fq = Fgy (m^,m2,m|,p^) and Fj = 
Fj'\mlmlmlp'), j = 1,2,3. 



3 Handling of the asymptotic expansion 

The asymptotic expansion at large p"^ values was proposed in [16], but only 
few terms were given there explicitly. It can be written symbolically as 

Fi = E,,o + log(p2)(E,,i) + log2(p2)(E,,2) , i = 0,1,2,3, (7) 

where Sjj- are power series in 1/p'f.. For the purpose of the program presented 
here, up to 18 terms in inverse powers of p^ were evaluated. That required 
the solution of a system of four linear equations per set of coefficients. Each 
coefficient in the expansion is a polynomial in three masses and logarithms 
of the masses. As the length of the coefficients is growing with the growing 
inverse power of p"^ we have made an effort to shorten the expressions using 
symmetric polynomials for Fq. The polynomials are calculated once and used 
several times in the evaluation of the coefficients. For the Fi {i = 1,2,3) we 
found out that the most effective way to simplify the coefficients is to make use 
of the relation Fi{n,m,1,m,l,m,l,p'^) = — (9Fo(n, mi,m2,m|,p^)/9mf after the 



simplification of tlie coefficients of Fq. All that resulted not only in shortening 
of the expressions, but what is more important, in significant (about 10 times) 
gain in CPU time necessary for the calculation of the sunrise master integrals. 

To speed up the convergence of the series and to enlarge the convergence region 
we use the Bernoulli change of variables, introduced in [20] and systematically 
used in [21,15], separately for each of the series Ejj from Eq.(7). Each of the 
series S (we drop here the subscript to shorten the expressions) 




S = 
after the (Bernoulli) change of variable 
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becomes 
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We obtain an expansion in y which, with 18 terms, provides with double 
precision (real*8) results for all values of p^ outside the interval [ -1.5,0.5] and 
for arbitrary masses. The precision of the result is estimated by taking the 
ratio of the last term to the sum of all the terms. As a matter of fact the 
double precision accuracy is obtained for particular values of masses also in a 
wider region of p^ (see next section). The estimation of the error was checked 
against the results of a Runge-Kutta method of comparable precision. It is to 
be recalled here that the convergence of the expansion in powers of 1/p^ is 
superior to the Runge-Kutta approach for large values of p^, so that the check 
is a really stringent one in the region —1.5 < pi < —1 and < p^ < 0.5 at the 
borderline of the convergence of the expansion. 

As the coefficients of the Bernoulli accelerated series are obtained numerically, 
we have checked the numerical stability of the procedure. The expected can- 
cellations occurring in Eq.(ll) were never affecting the accuracy of the result 



and the formulae remained numerically stable at variance with the formulae 
for the expansion at p^ = 0, which were thus not used in the program. 

For the values of p^ for which the asymptotic expansion does not work we 
use the direct numerical solution of the system of differential equations by 
means of the Runge-Kutta method described in details in [11]. The algorithm 
works relatively fast in the region of small pi values. Thus the combination 
of the methods, used in the present program, allows for the fast and accurate 
evaluation of the sunrise Master Integrals for all values of p^. 



4 The outline of the program 



The way the program works is shown schematically in Fig. 1. First, the pro- 
gram reads from the file input_BOKASUN.dat the values of p^, mi, m2, ma, the 
required relative accuracies (denoted as A^., A,. = (acc(l), acc(2), acc(3)) ) for 
the real parts, the imaginary parts and the moduli of the four functions, which 
the user wishes to calculate. The program checks if the value of the squared 
rescaled four momentum (p^) is within the interval A = [—1,0], where the 
series expansion is not valid. For p^. & A the program uses the Runge-Kutta 
method, described in [11], and gives the values of all functions -F„. The pro- 
gram tries to reach the accuracy asked by the user and gives the values with 
the best obtained accuracy, even if the required accuracy was not reached. 
Outside the interval A the program evaluates first the value of the functions 
Fn by using the asymptotic expansion described in the previous section. If 
the estimated relative accuracies for real and imaginary parts are lower than 
10^^^ or the accuracies are better (or equal) than (to) the accuracies asked by 
the user the program writes the F„ with their relative accuracies for real and 
imaginary parts. If any of the accuracy requirements is not met the program 
uses Runge-Kutta method to calculate all functions F„, n = 0,1,2,3. The 
result with better accuracy is written as the output. 

There are two subroutines to calculate the Mi's: 

bokasun, where the series expansion coefficients are calculated for each call, 
and 

bokasun_s, where the series expansion coefficients are calculated in advance 
in the subroutine prepare_store 

The second option is useful when one needs to calculate the Mi's for fixed 
masses and various p"^ values. It speeds up the calculations about 25 times. 

Both subroutines are called with parameters: 
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Fig. 1. The flowchart of the BOKASUN program (see the text for details). A^ (Ark)'- 
accuracy obtained by the expansion (Runge-Kutta) method. F^ {Fn^), n = 0, ...,3: 
results obtained by means of the expansion (Runge-Kutta) method. The A^ < A^fc 
condition is checked separately for the real and the imaginary part of each MI. Thus, 
when the Mi's are calculated by both methods, the result with better accuracy is 
given in the output. 

p2 , ml in , m2in , m3in , ace , Fn , delt are , deltaim 
declared as 



real*8 p2,mlin,m2in,m3in,acc(3) 

complex*16 Fn(0:3) 

real*8 deltare(0:3) ,deltaini(0:3) 



p2 is the square of the four momentum 

mlin,m2in,m3in are the internal masses 

acc(l) is the required relative precision for the real parts of the Mi's 

ace (2) is the required relative precision for the imaginary parts of the Mi's 

ace (3) is the required relative precision for the modulus of the Mi's; used 
only for Runge-Kutta method 

Fn(i) , 1=0,1,2,3 finite part of the ith MI 

deltare(l) ( deltaim(l)) relative accuracy of the real (imaginary) part of 
the ith MI 



5 Tests of the program and typical run times 



In [11] many tests of the Runge-Kutta method were performed in all regions 
of the p^ values. Comparisons have shown an excellent agreement between the 
code developed in [11] and the values published in [1,8]. In [12] the author 
states also the complete agreement with [11] of his code, published later as a 
part of [14]. In view of these comparisons we have just checked that the new 
part of the program, which uses the large pi expansion gives results which 
are in agreement with the Runge-Kutta method. As matter of fact, when 
the expansion in 1/p^ applies, our program reaches a precision of 10~^^ or 
better, which is higher than the other available programs, so that a direct 
comparison up to that precision was in general not possible. Nevertheless, 
extensive comparisons were made, limited to the relative precision of 10^^^ (or 
slightly better; that is the maximum precision which one can reach with the 
Runge-Kutta method for large pi values) for various sets of masses. The results 
were always in agreement within the errors. Machine precision comparisons 
were possible with [15] in the equal mass case, and an excellent agreement 
was found. 

The main gain in using the series expansion, whenever possible, is the reduc- 
tion of CPU time necessary to obtain the result. The CPU time (on a laptop 
with Intel Centrino Duo T7400 2.16 GHz processor) necessary for calculation 
of the Master Integrals with the double precision machine accuracy, with the 
expansion method, is about 8 ■ 10~^ s, which reduces to about 3 ■ 10^^ s when 
the expansion coefficients are calculated in advance. With the Runge-Kutta 
method it takes 10 s at p^ = 0.1 and 1800 s dX pi = 10 to obtain the relative 
accuracy of 10^^^. 



TEST RUN OUTPUT 



The distributed version of the program contains also a code of the test run, 
which uses both subroutines bokasun and bokasun_s. It provides also with 
an example of using the BOKASUN program. It reads the input parameters 
for the test run from the file input_BOKASUN.dat and appends the results 
of the Mi's with the obtained relative accuracies to files f O.dat (Fq), f 1 .dat 
(Fi), f2.dat (F2) and f3.dat (F3). In the files fO.dat.ref, fl.dat.ref, 
f2.dat.ref and f3.dat.ref, distributed together with the program, the re- 
sults expected for the test run are given. The test run takes about 6 minutes 
CPU as it calculates many points using the Runge-Kutta method. 
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